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A system of electrons in a local or nonlocal external potential can be studied with 1-matrix 
functional theory (1MFT), which is similar to density functional theory (DFT) but takes the one- 
particle reduced density matrix (1-matrix) instead of the density as its basic variable. Within 1MFT, 
Gilbert derived [PRB 12, 2111 (1975)] effective single-particle equations analogous to the Kohn- 
Sham (KS) equations in DFT. The self-consistent solution of these 1MFT-KS equations reproduces 
not only the density of the original electron system but also its 1-matrix. While in DFT it is 
usually possible to reproduce the density using KS orbitals with integer (0 or 1) occupancy, in 
1MFT reproducing the 1-matrix requires in general fractional occupancies. The variational principle 
implies that the KS eigenvalues of all fractionally occupied orbitals must collapse at self-consistency 
to a single level, equal to the chemical potential. We show that as a consequence of the degeneracy 
the iteration of the KS equations is intrinsically divergent. Fortunately, the level shifting method, 
commonly introduced in Hartree-Fock calculations, is always able to force convergence. We introduce 
an alternative derivation of the 1MFT-KS equations that allows control of the eigenvalue collapse by 
constraining the occupancies. As an explicit example, we apply the 1MFT-KS scheme to calculate 
the ground state 1-matrix of an exactly solvable two-site Hubbard model. 

PACS numbers: 71.15.Mb,31.15.xr 



I. INTRODUCTION 

Density functional theory (DFT) benefits from op- 
crating with the electron density, which as a function 
of just three coordinates is much easier to work with 
than the full many-body wavefunction. According to the 
Hohenberg-Kohn (HK) theorem^ the density of an elec- 
tron system in a local external potential v(r) may be 
found by minimizing a universal energy functional E v [n] , 
whose basic variable is the density. Remarkably, the den- 
sity uniquely determines the ground state wavefunction 
(if it is nondegenerate) , i.e., there can be only one ground 
state wavefunction yielding a given density, no matter 
what v{r) is. However, if the external potential is non- 
local, then the density alone is generally not sufficient to 
uniquely determine the ground state (see Appendix A for 
a simple example). GilbertP extended the HK theorem 
to systems with nonlocal and spin dependent external 
potential v(x,x'), where x — (r,a). It was proved that 
i) the ground state wavefunction is uniquely determined 
by the ground state 1-matrix (one-particle reduced den- 
sity matrix) and ii) there is a universal energy functional 
E v [7] of the 1-matrix, which attains its minimum at the 
ground state 1-matrix. The 1-matrix is defined as 

•y(x, x') = N J dx2 ■ ■ ■ dxNp(x, X2, ■ ■ ■ x^; x', X2, ■ ■ ■ xn), 

where J dx — J2a I anc ^ P = J2i w i O^il * s the 
full A^-clcctron density matrix with ensemble weights «',■ 
such that ^2;Wi = 1. An external potential may be 
nonlocal with respect to the space coordinates and/or 
the spin coordinates. For example, pseudopotentials are 
nonlocal in space, and Zeeman coupling — (h\e\ / mc) B ■ a , 



where a is the vector of Pauli matrices, is nonlocal in 
spin space. The coupling of electron motion and an ex- 
ternal vector potential, (\e\/2mc)(p ■ A + A ■ p), may also 
be treated as a nonlocal potential because p is a differen- 
tial operator. It is rather intuitive that for such external 
potentials, which couple to the system in more complex 
ways than the local potential v(f), it is necessary — in 
order to permit statements analogous to the HK theorem 
- to ref ine the basic variable accordingly. Hence spin- 
DFTj^l whose basic variables are the density and the 
magnetization density, applies to systems with Zeeman 
coupling. Current-DFTpSl whose basic variables are the 
density and the paramagnetic current density, has the 
scope to treat systems in which the current is coupled to 
an external magnetic field. Generally, if one considers an 
external potential that is nonlocal in space and spin, the 
necessary basic variable is the one-matrix^ which con- 
tains all of the single-particle information of the system, 
including the density, magnetization density and param- 
agnetic current density. 

The DFT-type approach that takes the 1-matrix as ba- 
sic variable will be referred to here as 1-matrix functional 
theory (1MFT). As in DFT, an exact and explicit energy 
functional is generally unknown. An important difference 
between 1MFT and DFT is that the kinetic energy is a 
simple linear functional of the 1-matrix, while it is not a 
known functional of the density. Thus, in 1MFT the only 
part of the energy not known explicitly is the electron- 
electron interaction energy W[7]. Several approximate 1- 
matrix energy functionals have been proposed and tested 
recently (see Refs. [7] and |5] and references therein.) No- 
tably, the so-called BBCrt approximations, 7 which are 
modifications of the Buijse-Baerends functional,^ have 
given fairly accurate results for the potential energy 
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curves of diatomic molecules^ and the momentum distri- 
bution and correlation energy of the homogeneous elec- 
tron gasP In Ref. [HI a density dependent fitting param- 
eter was introduced into the BBC1 functional such that 
the resulting functional yields the correct correlation en- 
ergy of the homogeneous electron gas at all values of den- 
sity. There is also the prospect of using 1MFT to obtain 
accurate estimates for the band gaps of non-highly corre- 
lated insulators!^! Many of the approximate functionals 
that have been proposed are similar to an early approxi- 
mation by MiillerEQ 

Actual calculations in 1MFT are more difficult than 
in DFT. The energy functional £^[7] must be minimized 
in a space of higher dimension because the 1-matrix is a 
more complex quantity than the density. In the calcu- 
lations cited above, the energy has been minimized di- 
rectly by standard methods, e.g., the conjugate gradient 
method. In DFT the energy is generally not minimized 
by such direct methods. Instead, the Kohn-Sham (KS) 
scheme^ provides an efficient way to find the ground 
state density. In this scheme, one introduces an aux- 
iliary system of N noninteracting electrons, called the 
KS system, which experiences an effective local potential 
v s (r). This effective potential is a functional of the den- 
sity such that the self-consistent ground state of the KS 
system reproduces the ground state density of the inter- 
acting system. It is interesting to ask whether there is 
also a KS scheme in 1MFT. The question may be stated 
as follows: does there exist a 1-matrix dependent effec- 
tive potential v s (x,x') such that, at self-consistency, a 
system of noninteracting electrons experiencing this po- 
tential reproduces the exact ground state 1-matrix of the 
interacting system? Although Gilbert derived such an 
effective potentialpl the implications were thought to be 
"paradoxical" because the KS system was found to have a 
high (probably infinite) degree of degeneracy. Evidently, 
the KS eigenvalues in 1MFT do not have the meaning of 
approximate single-particle energy levels, in contrast to 
DFT and other self-consistent-field theories, where the 
eigenvalues may often be interpreted as the negative of 
ionization energies, owing to Koopmans' theorem. The 
status of the KS s chem e in 1MFT appears to have re- 
mained unresolved) 13 ! 14 ! and recently it has bee n argued 
that the KS scheme does not exist in iMFTPiSIlIIH 
Gilbert derived the KS equations from the stationary 
principle for the energy. The KS potential was found 
to be 

v s (x, x') = v(x, x') + 6W/6i(x', x). (2) 

In this article, we propose an alternative derivation of the 
KS equations, which, in our view, gives insight into the 
nature of the "paradoxical" degeneracy of the KS system. 

One-matrix energy functionals are often expressed in 
terms of the so-called natural orbitals and occupation 
numbers.^ This makes them similar to "orbital depen- 
dent" functionals in DFT. The natural orbitals are the 
eigenfunctions of the 1-matrix, and the occupation num- 
bers are the corresponding eigenvalues.^ These quanti- 



ties play a central role in 1MFT. Recently, it was shown 
that when a given energy functional is expressed in terms 
of the natural orbitals and occupation numbers, the KS 
potential can be found by usin g a ch ain rule to evaluate 
the functional derivative in Eq. |2j 18 

Although the concept of the KS system can indeed 
be extended to 1MFT, it has in this setting some very 
unusual properties. In particular, the KS orbitals must 
be fractionally occupied, for otherwise the KS system 
could not reproduce the 1-matrix of the interacting sys- 
tem, which always has noninteger eigenvalues (occupa- 
tion numbers). This is different from the situation in 
DFT, where it is usually possible to reproduce the den- 
sity using only integer (0 or 1) occupation numbers, or 
in any case, only a finite number of fractionally occu- 
pied states. Due to the necessity of fractional occupation 
numbers, the 1MFT-KS system cannot be described by 
a single Slater determinant. However, we find that it can 
be described by an ensemble of Slater determinants, i.e., 
a mixed state. In order that the variational principle is 
not violated, all the states that comprise the ensemble 
must be degenerate. This implies that the eigenvalues of 
all fractionally occupied orbitals collapse to a single level, 
equal to the chemical potential. The degeneracy has im- 
portant consequences for the solution of the KS equations 
by iteration. We prove that the iteration of the KS equa- 
tions is intrinsically divergent because the KS system has 
a divergent response function Xs — Sj/Sv s at the ground 
state. Fortunately, convergence can always be obtained 
with the level shifting method.^ To illustrate explicitly 
the unique properties of the 1MFT-KS system, we apply 
it to a simple Hubbard model with two sites. The model 
describes approximately systems which have two local- 
ized orbitals with a strong on-site interaction, e.g., the 
hydrogen molecule with large internuclear separation.^ 
The Schrodingcr equation for this model is exactly solv- 
able, and we find that the KS equations in 1MFT and in 
DFT can be derived analytically. It is interesting to com- 
pare 1MFT and DFT in this context. We demonstrate 
that divergent behavior will appear also in DFT when 
the operator 1 — XsX -1 > where x an d Xs are the density 
response functions of the interacting and KS systems, re- 
spectively, has any eigenvalue with modulus greater than 
1 . In this expression the null space of x is assumed to be 
excluded. 

This article is organized as follows. In Sec.|TTJ we derive 
the KS equations in 1MFT and discuss how to solve them 
self-consistcntly by iteration. In Sec. |III| we compare 
three approaches to ground state quantum mechanics — 
direct solution of the Schrodinger equation, 1MFT and 
DFT — by using them to solve the two-site Hubbard 
model. 



II. KOHN-SHAM SYSTEM IN 1MFT 

It is not obvious that a KS-type scheme exists in 1MFT 
for the following reason. Recall that in DFT the KS sys- 
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tem consists of TV noninteracting particles and reproduces 
the density of the interacting system. The density of the 
KS system, if it is nondegenerate, is the sum of contri- 
butions of the TV lowest energy occupied orbitals 

occ 

n(f)=J2\<t>i(?)\ 2 - (3) 

i 

On the other hand, in 1MFT the KS system should repro- 
duce the 1-matrix of the interacting system. The eigen- 
functions of the 1-matrix are the so-called natural or- 
bitals, and the eigenvalues are the corresponding occupa- 
tion numbers.^ Occupying the TV lowest energy orbitals 
in analogy to pj, one obtains 

occ 

j(x,x')=Y,<Pi(x)<t>i(x')- (4) 

i 

Such an expression, in which the orbitals have only inte- 
ger (0 or 1) occupation, cannot reproduce the 1-matrix of 
an interacting system because the orbitals of an interact- 
ing system have generally fractional occupation (see the 
discussion in the following section.) The difference be- 
tween the 1-matrix in Q and the 1-matrix of an interact- 
ing system is clearly demonstrated by the so-called idem- 
potency property. The 1-matrix in Q is idempotent, i.e., 
J dx ,r y(x, x') , y(x', x") — 7(2,2;"), while the 1-matrix of 
an interacting system is never idempotent. However, if 
the KS system is degenerate and its ground state is an 
ensemble state, the 1-matrix becomes 

~f(x,x') = J2fi<f>i(x)<f>*i(x') (5) 

i 

with fractional occupation numbers The TV-particle 
ground state density matrix of the KS system is p s — 
Yli w i l^i) where the <&i are Slater determinants 

each formed from TV degenerate KS orbitals. The oc- 
cupation numbers fi are related to the ensemble weights 
m by 

fi = ^2 w j®ji, ( 6 ) 
3 

where Qji equals 1 if <\>i is one of the orbitals in the 
determinant <i> 7 - and otherwise. 21 



A. Derivation of the 1MFT Kohn-Sham equations 

In this section, we discuss Gilbert's derivatiorP of the 
KS equations in 1MFT and propose an alternative deriva- 
tion. We begin by reviewing the definition of the univer- 
sal 1-matrix energy functional £^[7]. 

One-matrix functional theory describes the ground 
state of a system of TV electrons with the Hamiltonian 

H = Y!Lx{U + «i) + W, where * = ~V 2 /2 is the kinetic 
energy operator, v is the local or nonlocal external poten- 
tial operator, and W = J2i<j fai ~ is the electron- 
electron interaction (in atomic units h = m = e = 1). 



The ground state 1-matrix and ground state energy can 
be found by minimizing the functional 

E v [ 1 ]=Tr((t + v) 7 ) + W[~ f ], (7) 

where 

W[i\ = {^o\W\^o)- (8) 

By extending the HK theorem, Gilbert provecP that 
a nondegenerate ground state wavefunction, ^> Q , is 
uniquely determined by the ground state 1-matrix, i.e., 
^0 is a functional of 7. For this reason the interaction en- 
ergy, as defined in Q, is a functional of 7. It is apparent 
that ^ defines W[j] only for 7 that are the ground state 
1-matrices of some system (with Hamiltonian H). In 
this article, a 1-matrix is said to be w-representable (VR) 
if it is the ground state 1-matrix of some system with 
local or nonlocal external potential. Gilbert remarked 
(see the discussion between equations (2.24) and (2.25) 
in Ref. |2j that, in principle, the domain of W[j] can 
be extended to the space of ensemble TV-representablc 
(ENR) 1-matrices. A 1-matrix is said to be ENR if it 
can be constructed via (fTl) from some TV-particle den- 
sity matrix p = Wi which is not required to 
be a ground state ensemble. One possible extension to 
the ENR space is pro vided by the so-called constrained 
search functionaP^U 

W[j] = mm^Tr(Wp), (9) 

where the interaction energy Tr(Wp) is minimized in the 
space of TV-particle density matrices p that yield 7 via 
([lj . The definition ^ is a natural extension to the ENR 
space because when it is adopted ^ may be expressed 
as 

E v [-y]=mm^Tr{Hp). (10) 

This is a variational functional which attains its mini- 
mum at the ground state 1-matrix, as seen from 

min 7 £J„[7] = mm^Tr(Hp) — E , (11) 

where Eq is the ground state energy. The extension to 
the ENR domain is significant, especially for applications 
of the variational principle, because the conditions a 1- 
matrix must satisfy to be ENR are known and simple to 
impose on a trial 1-matrix, while the conditions for v- 
representability are unknown in general. The necessary 
and sufficient conditions^ for a 1-matrix 7 to be ENR 
are i) 7 must be Hermitean, ii) J dx~f(x, x) = TV, and iii) 
all eigenvalues of 7 (occupation numbers) must lie in the 
interval [0,1]. The third condition is a consequence of 
the Pauli exclusion principle. 

The 1MFT-KS equations were derivecP from the sta- 
tionary conditions for the energy with respect to arbi- 
trary independent variations of the natural orbitals 4>i 
and angle variables 6i chosen to parametrize the occupa- 
tion numbers according to fi = cos 2 9i (0 < 0i < tt/2). 
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For the purpose of describing a variation in the ENR 
space, this set of variables, namely {5<fii,8(fr*. l 56i}, is re- 
dundant. An arbitrary set of such variations may or may 
not correspond to an ENR variation, and when it does the 
variations will not be linearly independent. This causes 
no difficulty, of course, because it is always possible to for- 
mulate stationary conditions in a space whose dimension 
is higher than necessary, provided the appropriate con- 
straints are enforced with Lagrange multipliers. Accord- 
ingly, the Lagrange multiplier terms J2ij ^ij{{<t>i\4>j) ~^ij) 
which maintain the orthogonality of the orbitals and the 
term ptjY], fi — N) which maintains the total particle 
number were introduced. The KS equations were found 
to be 



(i+v s ) \<j>i 



(12) 



where the kernel of the effective potential is v s {x,x') — 
v(x,x') + 8W/S'y(x' ,x), if the functional derivative ex- 
ists. The stationary conditions imply that all fractionally 
occupied KS orbitals have the same eigenvalue e, = /i. 
Gilbert described this result as "paradoxical" because in 
interacting systems essentially all orbitals are fractionally 
occupied. 

The above stationary conditions assume E v to be sta- 
tionary with respect to variations in the ENR space (ex- 
cept variations of occupation numbers equal to exactly 
or 1 in the ground state, which are excluded by the 
parametrization) . However, it is not known, in general, 
whether E v is stationary in the ENR space; the mini- 
mum property (II]) ensures only that it is variational. 
Recall that the ENR space consists of all 7 that can be 
constructed from an ensemble, and the energy of an en- 
semble is not stationary with respect to variations of the 
many-body density matrix p. Therefore, the stationary 
conditions applied in the ENR space may be too strong. 
On the other hand, the quantum mechanical variational 
principle guarantees that E v is stationary with respect 
to variations in the VR spaced but it is not known how 
to determine whether a given 7 is VR. Hence, it is not 
known how to constrain the variations of 7 to the VR 
space. Nevertheless, it may be that in some systems the 
entire neighborhood of j gs in the ENR space is also VR. 
In such cases, E v is stationary in the ENR space, and 
the stationary conditions applied in Ref. 2 are satisfied 
at the ground state. 

We find it helpful to construct an alternative derivation 
of the KS equations. Consider the energy functional 



Qvh] = Kb] - e i(f} ~ 



(13) 



where ej — £j({qk}) are Lagrange multipliers that con- 
strain the occupation numbers fj of the natural orbitals 
to chosen values qj, which satisfy < qj < 1 and 
J2j qj = N. These Lagrange multipliers allow us to in- 
vestigate the degeneracy of the KS eigenvalues, which 
leads to the "paradox" described by Gilbert. We have 
omitted the Lagrange multipliers Ay used in Gilbert's 



derivation. They are not necessary in our derivation be- 
cause we formulate the stationary conditions with respect 
to variations of the 1-matrix instead of the orbitals. We 
adopt the definition (JsJ) for IF [7], so the domain of Gv[l] 
is the ENR space. A variation £7 will be said to be ad- 
missible if j gs + 8j is ENR. For convenience we assume 
that the static response function x — 8^/dv for the in- 
teracting system under consideration has no null vectors 
apart from the null vectors associated with a) a constant 
shift of the potential (which is a null vector also in DFT) 
and b) integer occupied orbitals, i.e., orbitals with occu- 
pation numbers exactly or 1. If \ has additional null 
vectors, the following derivation must be modified; the 
necessary modifications are discussed below. Granting 
the above assumption, Q v is guaranteed to be stationary, 
and the KS equations can be derived from the stationary 
condition 8Q V = with respect to an arbitrary admissible 
variation of 7. The first variation of Q v is 



8Q V 



Tr((i + v)5j) + Tr(wSj) - ^ e 3 Sf 3 

j 

^2(<f>i \ (i + v + w)\ <t>j) (^l^l^i) 

ij 

-J2 e j {<i>Mi\<i>i) 

ij 

^2(h i:j - ejSij) Sjji, (14) 



where the variation of the 1-matrix is expressed as 8"/ij = 
(<j>i #7 in the basis of the ground state natural or- 
bitals, and the relation 8W = Tr(w8j) defines a single- 
particle operator w. In ( 14), we have also introduced the 



definition of the Hermitean operator 



t + v + w, 



(15) 



which will be seen to be the KS Hamiltonian. If the 



last line of ( 14 1 is to be zero for an arbitrary Hcrmitian 



CjSij = for all i 



matrix Sj, then we must have hij 
and j. ENR condition (i) has been maintained explicitly 
by requiring the variation to be Hcrmitian. ENR condi- 
tions (ii) and (iii) do not impose any constraint on the 
space of admissible variations as they are maintained by 
the Lagrange multipliers e^. The matrix elements are 
functionals of the 1-matrix, and the 1-matrix that satis- 
fies the stationary conditions — £j5ij = can be found 
by solving self-consistently the single-particle equations 



h\<l>i) 



i) 



(16) 



together with ([5]). These are the KS equations in 1MFT. 
If they are solved self-consistently with the occupation 
numbers fixed to the values qi, they give the orbitals 
which minimize E v subject to the constraints /, = qi. 
The KS potential is v s = v + w. The term w is the 
effective contribution of the electron-electron interaction 
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to the KS potential. In coordinate space, its kernel is 



){x,x') — ( x I w\x') 
SW 
S^f '{x 1 ', x) ' 



(17) 



which recovers Gilbert's result. The kernel of the KS 
Hamiltonian may be written in the familiar form 

h(x,x') = 5(x-x')(-^l)+v(x,x') 

+5(x - x')vh(x) + v xc (x,x'), (18) 

where v(x, x') is the external potential and w(x,x') 
has been divided into the Hartree vh(x) and exchange- 
correlation v xc (x,x') potentials. In 1MFT, the exchange- 
correlation potential is nonlocal. 

The 1MFT-KS scheme optimizes the orbitals for a cho- 
sen set of occupation numbers, but it does not itself pro- 
vide a rule for choosing the occupation numbers. On this 
point, it is different from the DFT-KS scheme, where 
the occupation numbers are usually uniquely determined 
by the aufbau principle (T = Fermi statistics).^ In 
1MFT, the KS equations have a self-consistent solution 
for any chosen set of occupation numbers {qi} that sat- 
isfy < qi < 1 and J2i 1i = N. The Lagrange multipliers 
Cj, which are seen to be the KS eigenvalues, adopt val- 
ues such that the minimum of Q v occurs for a 1-matrix 
Imin whose occupation numbers are precisely the set 
{qi}. Therefore, the unconstrained minimum of Q v coin- 
cides with the minimum of E v subject to the constraints 
fi = Qi- To find the ground state occupation numbers, 
{qf s }, one can search for the minimum of the function 
G v ({qi}) = min 7 £/„[7]> where the minimization of Q v [j] 
can be performed by the KS scheme. 

What can be said about the KS eigenvalues e^? As the 
minimum of Q v is a stationary point, we have 







dQv_ 
dfj 



dE v 



ejitii}) (19) 



for all j. E v is a variational functional in the ENR space. 
It attains its minimum at the ground state 1-matrix "f gs , 
where dE v /dfj must vanish for all fractionally occupied 
(0 < fj < 1) orbitals, for otherwise the energy could be 
lowered. When = q^ s for all i, ^ m i n = 7 gs , and ( 19 1 im- 
plies ej({gf s }) = for all fractionally occupied orbitals. 
Thus, we find that the KS eigenvalues of all orbitals that 
are fractionally occupied in the ground state must col- 
lapse to a single level when the chosen set of occupa- 
tion numbers approach their ground state values, i.e., as 
qi — > qf s . In Gilbert's derivation of the 1MFT-KS equa- 
tions, all fractionally occupied KS orbitals were found to 
have the eigenvalue gj — /i, where /1 is the chemical po- 
tential. As we have not introduced the chemical potential 
in our derivation (we consider a system with a fixed num- 
ber of electrons) , the eigenvalues collapse to instead of 
fi. The above arguments do not apply to orbitals with 
occupation numbers exactly or 1 because these values 



lie on the boundary of the allowed interval [0, 1] specified 
by ENR condition (iii). All that can be concluded from 
the fact that E v has a minimum in the ENR space is 
€j > for orbitals with fj = and ej < for orbitals 
with fj = 1. States with occupatio n num bers exactly or 
1 have been called "pinned states. "fHUHI Instances of such 
states in real systems have been reported f23 though their 
occurrence is generally considered to be exceptional ISliflJ 
Due to the collapse of the eigenvalues at the ground 
state, the KS Hamiltonian becomes the null operator 



(20) 



in the subspace of fractionally occupied orbitals. This is 
of course analogous to the familiar condition df /dx = 
for the extremum of a function f(x). Gilbert describecP 
a similar result (with fj,I replacin g 0) a s "paradoxical," a 
statement that has been repeated ! 13 l 14 l The problem with 



( 20 1 is that while we expect the KS Hamiltonian to de- 



fine the natural orbitals, any state is an eigenstate of the 
null operator. However, the KS Hamiltonian is a func- 
tional of the 1-matrix, and when the occupation num- 
bers are perturbed from their ground state values, the 
degeneracy is lifted and the KS Hamiltonian does define 
unique orbitals. In the KS scheme outlined above, this 
corresponds to the optimization of the orbitals with oc- 
cupation numbers fixed to values qi, perturbed from the 
ground states values. In the limit that the occupation 
numbers approach their ground state values, the opti- 
mal orbitals approach the ground state natural orbitals. 
The degenerate eigenvalues generally split linearly with 
respect to perturbations away from the ground state. In 
particular, 



dei 
dqj 



= / dydy' 1 



5h d-f(y,y') 



10.) 



8j(y, y') dq.j 
= - J dxdx'dydy'(j)*(x)(t)*(y')x^ 1 (xx',yy') 
x^(x')cf>j(y). (21) 



Here, x is the static response function defined as 



X{x,x';y,y') = 



S r y(x, x 1 ) 
Sv(y,y') ' 



(22) 



The relation Sh/Sj = — x~ X used in (21) is derived in 
the following section. If \ has a null space, its inverse 
is defined only on a restricted space. For example, (21 1 



does not apply to pinned states as there is a null vector 
associated with each pinned state (see below). 

Our derivation of the 1MFT-KS equations in fact as- 
sumes that the static response function \ of the interact- 
ing system has no null vectors except for those associated 
with pinned states and a constant shift of the potential. 
We now show that this guarantees Q v is stationary. If 
the interacting system has any other null vectors we can 
no longer be certain that Q v is stationary and the deriva- 
tion should be modified as described below. We have 
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remarked already (see Ref. 34) that Q v is stationary 
in the VR space, i.e., it satisfies the stationary condi- 
tion SGv — with respect to an arbitrary variation of 
the 1-matrix in the VR space. However, our derivation 
of the KS equations requires Q v to be stationary in the 
ENR space. As the VR space is a subspace of the ENR 
space, this is a stronger condition. The assumption that 
X has no null vectors (apart from those associated with 
pinned states) is equivalent to assuming that any ENR 
variation (apart from variations of the pinned occupa- 
tion numbers) is also a VR variation. For if \ has no 
null vectors, then it is invertible and any ENR variation 
<$7 can be induced by the perturbation Sv = x -1 ^.^ 
Hence, with the above assumption, Q v is guaranteed to 
be stationary with respect to any ENR variation. The 
above arguments do not apply to variations of pinned 
occupation numbers because there are null vectors asso- 
ciated with such variations; nevertheless, Q v is station- 
ary with respect to such variations as this is maintained 
by the Lagrange multipliers e^. It is now clear how to 
modify the derivation of the KS equations when \ has 
additional null vectors. By introducing new Lagrange 
multipliers, the additional null vectors can be treated in 
analogy with the pinned states. For example, suppose 
X has one additional null vector u = J2ij u ij \4>i) {4>j\-< 
where u is hermitian and Tr{uu) = 1. The new energy 
functional Q' v [y\ — Q v \l\ ~ k{1u — Pu) will be station- 
ary with respect to an arbitrary variation in the ENR 
space. Here, the Lagrange multiplier k enforces the con- 
straint 7,j = p Ul where 7„ = Tr{^u) is the component of 
7 corresponding to it. The stationary condition SQ' V — 



leads to the set of equations hi 



in 



the basis of natural orbitals. The 1-matrix that satisfies 
these equations can be found by solving self-consistently 
the eigenvalue equation h = u>i together with 

j{x,x') = T,ijk1i^ji S ki^j( x )^k{ x ')> where s is the uni- 
tary matrix that diagonalizes the matrix u in the basis of 
natural orbitals, i.e., SuS^ is diagonal. The energy of the 
self-consistent solution defines a function G' v ({qi} , p u ) , 
whose minimum with respect to {%} and p u is the ground 
state energy. It may not be known in advance whether 
the response function of a given interacting system will 
have null vectors. Therefore, it is helpful to understand 
how null vectors occur. 

Null vectors of x a re conn ected with the so-called 
nonuniqueness problerri- 3 ! 25 ! 26 ! in various extensions of 
DFT. A system with the ground state is said to have 
a nonuniqueness problem if there is more than one ex- 
ternal potential for which ^ is the ground state. The 
Schrodinger equation defines a unique map from the ex- 
ternal potential to the ground state wavefunction (if it is 
nondegenerate) , but when there is more than one exter- 
nal potential yielding the same ground state wavefunc- 
tion, the map cannot be inverted. In 1MFT the general- 
ity of the external potential (nonlocal in space and spin 
coordinates) allows greater scope for nonuniqueness than 
in the other extensions of DFT. Of course, every degree 
of nonuniqueness is a null vector of x because if Sv does 



not change the ground state wavefunction, it does not 
change the 1-matrix either, and hence it is a null vector. 
In fact, every null vector of x is caused by nonuniqueness; 
the existence of a null vector Sv that induced a nonzero 
S^o would contradict the one-to-one relationship 7 <-> ^0 
proved by the extension of the HK theorem to 1MFTP 
It was mentioned above that there are null vectors asso- 
ciated with the pinned states. Suppose <pk is a natural 
orbital with occupation number fk = in the ground 
state. The perturbation of the external potential Sv = 
A |</>/.) (</>fc| does not change the ground state if the system 
has an energy gap between the ground state and excited 
states and A is small enough because SV't'o = 0, where 
SV = J dxdx'^(x)Sv(x,x')tp(x') and ip and ^ are field 

operators. If fk = 1, SV^fo = ^0 an d the ground state 
is again unchanged by the perturbation. The "vector" 
\4>k) (0fc| is therefore a null vector of x if ^fe is a pinned 
state. Another type of nonuniqueness, which has been 
called systematic nonuniqueness J^Sl is related to constants 
of the motion. Suppose A = J dxdx'tjj 1 ' (x)a(x, x')ip(x') is 
a constant of the motion. The ground state, if it is nonde- 
generate, is an eigenstate of A as constants of the motion 
commute with the Hamiltonian. If the system has an en- 
ergy gap between the ground state and the first excited 
state, then a perturbation SV = XA will not change the 
ground state wavefunction if A is small enough. Thus, A 
is a null vector of x- 

As in DFT, an exact and explicit expression for the 
universal energy functional in 1MFT is unknown in gen- 
eral. In actual calculations it is usually necessary to use 
approximate functionals. Many of the approximate en- 
ergy functionals that have been introduced are expressed 
in terms of the natural orbitals and occupation numbers. 
Such functionals are valid 1-matrix energy functionals, 
but as the dependence on the 1-matrix is implicit rather 
than explicit, they have been called "implicit" function- 
als. Recently, the KS equations were derived for this 
case.^ It was found that the contribution to the KS po- 
tential from electron-electron interactions can be evalu- 



ated by applying the following chain rule to (17) 



w(x, x) 



SW 



5"f(x' , x) 



SW 6<j>i{y) 
Scj)i(y) Sj(x',x) 

SW S<f>*(y) 

H*{y) Sj(x',x) 

dW Sfj 
dfi Sj(x',x)' 



i J 

i 

E 



(23) 



B. Iteration of the KS equations 

In this section we show that the "straightforward" pro- 
cedure for iterating the KS equations ([5| and (16 1 is in- 
trinsically divergent. 
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The KS equations are nonlinear because the KS Hamil- 
tonian itself depends on the 1-matrix. In favorable cases 
such nonlinear equations can be solved by iteration. 
Given a good initial guess for the 1-matrix, iteration may 
lead to the self-consistent solution corresponding to the 
ground state. In order to iterate (16), one needs an al- 
gorithm to define the 1-matrix of iteration step n + 1 
from the 1-matrix of step n, i.e., one needs to "close" 
the KS equations. In the previous section, we saw that 
in the 1MFT-KS scheme the occupation numbers fi are 
held fixed during the optimization of the natural orbitals. 
The following is a "straightforward" algorithm that op- 
timizes the natural orbitals: i) the KS Hamiltonian for 
step n + 1 is defined by 



(n)l 



(24) 



where h is given by (151 and 7'") is the 1-matrix of it- 



eration step n; ii) the eigenstates of hy l+1 > are taken as 
the natural orbitals of step n+ 1; iii) the 1-matrix of step 
n + 1 is constructed from the natural orbitals of step n+ 1 
by the expression 
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(n+l) 



0, x') = y^fi 



(x)^ n+1 \x f ). (25) 



Let ui be the eigenstates of the KS Hamiltonian h^ n+1 \ 
In operation (ii), the natural orbitals 0[ n+1 ^ are chosen 
from among the Ui such that has maximum over- 

lap with (f)\ n \ i.e., Tr((^ n+1 ^> — j^) 2 ) is the minimum 
possible. If this procedure converges to the stationary 
1-matrix 7 m i„ giving the lowest energy possible for the 
fixed set of occupation numbers, then it defines the func- 
tion G v ({fi}), introduced in the previous section, for 
which the ground state energy is the absolute minimum. 
Unfortunately, for any set of occupation numbers {fi} 
sufficiently close to the ground state occupation num- 
bers, this procedure does not converge to j m in- In other 
words, the "straightforward" algorithm defines an itera- 
tion map for which the ground state 7 gs is an unstable 
fixed point. This is a consequence of the degeneracy of 
the KS spectrum at the ground state. 

The divergence of the iteration map is revealed by a 
linear analysis of the fixed point. Suppose the occupation 
numbers are fixed to values perturbed from their ground 
state values by Sfi. Let us consider an iteration step n 
and ask whether the next iteration takes us closer to the 
stationary point 7 m i„ that gives the minimum energy for 
the fixed occupation numbers. The linearization of the 
iteration map at the stationary point gives 



(n+l) 



;,(»+!) 



Xshmin](vi n+1) -V? in ) 



(26) 



Xs(x,x';y,y') 



6j(x, x') 

fi — fj 



EE 



x ^(x)ti(x'WAy)My')- (27) 



In the last line of ( 26 ) , we have used 



- 1 ) - h* 



-x #7 



(n) 



(28) 



which can be established by the following arguments. 
First consider 



i(n+l) _ jjinin 



6h[j\ 



57 



(<>) 



(29) 



which follows from (24 1. The KS Hamiltonian is an im- 



plicit functional of the 1-matrix, and (29) defines the 



first order change of the KS Hamiltonian with respect 
to a perturbation of that 1-matrix. Since the occupation 
numbers are close to their ground state values, we may 
make the replacement 



6h 

57 



5h 



(30) 



7gs 



which is valid to 0(max(|<5/j|)). Thus, to establish ( 28 1 



we need to show 8h/5~f — — x at the ground state 1- 
matrix ~/ gs . The KS Hamiltonian h is associated with 
the original many-body Hamiltonian H , which has ex- 



ternal potential v(x, x'). According to (20 P3, h[jg S ] = 0. 
Consider now a Hamiltonian H' with a slightly differ- 
ent external potential v'(x, x') — v{x, x') + 8v(x, x') such 
that its ground state 1-matrix is ~f' gs = j gs + #7. The 

associated KS Hamiltonian is h' — h + v' — v. At the new 
ground state, ^'[7^] — 0. This allows us to relate 8h to 
Sv as 



Sh = h[% s ] - h[ 7gs ] 

= Hlgs] + - M 



-Sv. 



(31) 



Finally, using Sv — 
which verifies ( 28 ) . 



X 1 Sj we obtain Sh/Sj = — x 



Returning to the question of convergence, we see that 
( 26 1 implies that the next iteration takes us farther 
from the stationary point 'Jrnin- The reason is that 
|detx s x _1 | > 1 if 7 is sufficiently close to the ground 
state. According to (211 the KS response diverges as 
7 - * Igs because — Cj ~ 0(max(|<5/;|)). For a fixed set 
of occupation numbers sufficiently close to their ground 
state values, the moduli of all eigenvalues of the opera- 



where v s = v+vh+v xc is the KS potential. The response tor XsX become greater than 1 (the null space of x 1S 
function \ was defined in (22 1. The KS response function assumed to be excluded). Therefore, a perturbation £7 
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from the ground state is amplified by iteration, and the 
ground state is an unstable fixed point of the iteration 
map. A fixed point is stable if and only if all eigenvalues 
of the linearized iteration map have modulus less than 1. 



C. Level shifting method 

In the previous section, it was shown that the 
"straightforward" iteration of the KS equations is intrin- 
sically divergent. To obtain a practical KS scheme the 
iteration map must be modified. In this section, we con- 
sider the level shifting method^ and by linearizing the 
modified iteration map we obtain a criterion for conver- 
gence. 

Intrinsic divergent behavior can be encountered also 
in the Hartree-Fock approximation, and various modifi- 
cations of the iteration procedure have been introduced, 
for example, Hartree damping (also called configuration 
mixing) and level shifting. The level shifting method is 
particularly attractive in 1MFT because it can prevent 
the collapse of eigenvalues that is the origin of divergent 
behavior. Indeed, it has been shown that the level shift- 
ing method is capable of giving a convergent KS scheme 
in 1MFTP 

In the "straightforward" iteration procedure, the 
change of the orbitals, to first order, from iteration step 
n to iteration step n + 1 is 



From the last line of (33 1, we obtain a criterion for the 



A 



2^ — t. ^'W' ( 32 ) 



where M™-* is the KS Hamiltonian for iteration step n 
and in the second line the orbitals and eigenvalues are 
from iteration step n. In the level shifting method, the 



first order change in the orbitals given by (32 1 is altered 
by applying the shifts — > e* + Q to the eigenvalues 
in the denominator. To first order, this modification is 
equivalent to adding the term A = J2i Cil^^i&i \ to 
the KS Hamiltonian for step n + 1. Let = h + A 
define the level shifted Hamiltonian. Repeating the linear 
analysis of the previous section for the iteration map for 
this level shifted Hamiltonian, we find 



6", 



(n+l) 



(n) 



(33) 



where we have defined the operator Cl with the kernel 
fl(xx\ yy') = f dzdz' Xs(xx' , zz') 



him/) 



convergence of the iteration map. All of the eigenvalues 
of the operator 



A = -X s [lmin]X 



(35) 



must have modulus less than 1. The dependence on the 
level shift parameters Q enters only through the shifted 
eigenvalues in the denominator of \s- The level shifting 
method is effective because it prevents the divergence of 
Xs at the ground state and there is a cancellation between 
the two terms in ( 35 1 . Unfortunately, the convergence 



criterion depends on x, which is unknown at the outset 
of a 1MFT calculation. In Sec. |TTT] the level shifting 
method is applied in an explicit example and the above 
criterion is verified. 



D. Properties of the KS system 

The distinguishing feature of the KS system in 1MFT 
is the degeneracy of the eigenvalue spectrum. This has 
surprising consequences. It was shown in section |II A| 
that the KS eigenvalue spectrum splits linearly as we 
move away from the ground state 1-matrix. There- 
fore, the total KS energy changes linearly with respect 
to the displacement, i.e., E s [^] — E s ["/ gs ] cx 5j, where 
£^[7] = tr(h[ j]'f) (for a specific example see Fig. [4] 



in Sec. 1KB 2 1. This is surprising because such linear 



changes do not occur for the energy functional E v (in 
the VR space) . The immediate implication is that E s [7] 
is not stationary at the ground state. While this causes 
no difficultly in principle — we need only the functional 
E v to be stationary — it is intimately connected with the 
divergence of the iteration map. Precisely at the ground 
state Esljgs] = X^ e «i where the prime indicates that 
only the pinned states with /i = 1 contribute to the sum. 
Away from the ground state the KS eigenvalue spectrum 
splits, and £^[7] is a multivalued functional due to the 
choice implied in occupying the new KS levels. This is the 
same choice e ncoun tered in the iteration of the KS equa- 
tions (see Sec. 



II B I, where the natural orbitals 



are 



selected from among the eigenstates of the KS Hamilto- 
nian. Near the self-consistent solution, there will be one 
such choice for which the resulting ^ n+1 '> is very close to 

It was shown in the preceding section that the static re- 
sponse function of the KS system diverges at the ground 
state. Thus, even an infinitesimal perturbation Sv s may 
induce a finite change of 7. At the ground state, all of 
the natural orbitals, except those which have an occu- 
pation number that is degenerate, are uniquely defined. 
The natural orbitals which belong to a degenerate occu- 
pation number are only defined modulo unitary rotation 
in the degenerate subspace. When a perturbation is in- 
troduced, the natural orbitals change discontinuously to 
the eigenstates of the perturbed KS Hamiltonian h = 8v. 
These eigenstates may be any functions in the degenerate 
Hilbert space because Sv is arbitrary. 
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III. 



TWO-SITE HUBBARD MODEL 



The 1MFT-KS system has some unusual features, such 
as the collapse of the KS eigenvalues at the ground state, 
so it is desirable to derive explicitly the KS equations for a 
simple model. The Hubbard model on two sites provides 
a convenient example because it is exactly solvable and 
especially easy to interpret. Also, analytic expressions for 
the 1-matrix energy functional and KS Hamiltonian can 
be obtained. In the following sections, for the purpose 
of comparison, we find the ground state of the two-site 
Hubbard model by three methods — direct solution of 
the Schrodinger equation, 1MFT and DFT. 



Direct solution 



in the basis ($0, $2, $3)- We have defined the dimension- 
less variable v = V/B. The secular equation \K — «,-J = 
is 



«f - (x 2 + l)n 2 + (x 2 - v 2 )m + v 2 y 2 = 0. 
The normalized eigenvectors of H are 




(39) 



(40) 



and have energy Ei — A + Bn^ for i = 0, 2, 3, where «j is 
a root of the secular equation ( 39 1 . We have also defined 



Pi = Ki(K» - x 2 ) ~ v 2 y 2 



(41) 



and 



The Hamiltonian of the two-site Hubbard model is 
H = f + U + V with 



U = U (hifhii + nafna;) 
V = v\ («! - n 2 ) , 



where t\% = <2i = t, cj CT and Ci„ are the creation and 
annihilation operators of an electron at site i with spin 
<7, and hi = Yl<T C i<r c icr- We consider only the sector of 
states with N — 2 and S z — 0, i.e., a spin unpolarized 
system. In this sector, the eigenstates of T + U are 





(36) 



$0 



$ 2 = 




(37) 



in the site basis {cj 



''IT 2i 



10), 



4 c t 

-2t 1-1 



|0). 



}. The following variables have been intro- 
duced x = cos(7r/4 — a /2), y = sin(-7r/4 — do/2), and 
tanao = U/4t with < ao < ir/2. The eigenvalues of 
T + U for the states are 



Ao 



By 2 , 



A 2 - B(x 2 - y 2 



Ai = 
A 3 





-Bx 2 



where B = \/U 2 + T 2 and T = At. $0, $2, and $3 are 
singlet states (S = 0) and $1 is a triplet state with S = 1 
and 5 2 = 0. We will omit $1 from consideration as it is 
not coupled to the other states by the spin-independent 
external potential chosen in ( 36 1 . The Hamiltonian may 
be written as H = Xnl + BK, where 



K 




(38) 



m = [ni{v 2 {Zx 2 -l) + y 2 ) + K l {x 2 y z {2v 2 - I)) 

,2„.2/,,2 „.2^ 1 / 2 



1 2 2/ 2 2\1 - 

+1/ y (y -y)\ 



(42) 



The two dimensionless energy scales of the system are 
the interaction strength U /T and the bias V/T. The 
behavior of the system with respect to these energy scales 
is illustrated in Fig. [l] The quantity m = (rii — n 2 )/2, 
where rn is the average ground state occupancy of site 
i, is plotted with respect to the external potential V for 
various values of the interaction strength U. For the 




FIG. 1: [color online] The density variable m = (nj — n-i) /2, 
is shown with respect to the dimensionless external poten- 
tial v = V/B (B = VT 2 + U 2 ). Curves for U/T = 
(1/16,1/4,1,4) are shown as (solid [blue], dotted [green], 
dash-dotted [light blue], dashed [red]) curves, respectively. 



ground state, 



(* H *o> 

2i>k 2 i x 2 , 



Vo 



K - X 



(43) 
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where rh — [h\ — n<i)/2. A weakly interacting system 
(e.g., the solid [blue] curve in Fig. [I]) responds strongly to 
the external potential. In contrast, a strongly interacting 
system (e.g., the dashed [red] curve) responds weakly up 
to a threshold V/B ~ 1 (for a strongly interacting system 
B ~ U.) This behavior has a simple interpretation: in 
order for the external bias to induce charge transfer, it 
must overcome the on-site Hubbard interaction. In the 
limit U — » oo, the curve develops step-like behavior near 
V/B ~ ±1. 



B. Solution by 1MFT 

In the first part of this section, we derive the en- 
ergy functional and KS Hamiltonian. In the second, 
we demonstrate the divergence of the iteration of the 
KS equations. In the third, we use the level shifting 
method^ to obtain a convergent KS scheme. 



1. Energy functional and KS Hamiltonian 

For lattice models such as the Hubbard model, the 1- 
matrix is defined as 



7(?ct,jt) = (#|4c jV |*) 



(44) 



One may ask whether the HK theorem (or Gilbert's ex- 
tension in 1MFT) applies when the density (or 1-matrix) 
is defined over a discrete set of points, i.e., when the 
continuous density function n(r) is replaced by th e site 
occupation numbers n,. This has been investigated) 27 1 28 1 
and it was found that the HK theorem remains valid. 
We consider here only spin unpolarized states (S z = 0). 
Accordingly, we define the spatial 1-matrix 

= ^2,7(i<r,3<r)- (45) 

a 

The 1-matrix may be expressed as, cf. ([5]), 

7(ij)=X)/«0aW&(j), (46) 



where 4> a are the spatial natural orbitals. As our sys- 
tem is spin unpolarized, the spin up and spin down spin- 
orbitals have the same spatial factors. Therefore, in (46) 
each spatial orbital 



b a may be occupied twice (once by a 
spin up electron and once by a spin down electron), i.e., 
< f a < 2. It is convenient to parametrize the natural 
orbitals as 



cos(0/2) 
sin(0/2) 



and 



sin(0/2) 
-cos(0/2) 



(47) 



In terms of this parametrization, the 1-matrix in the site 
basis is 



7 = I + A (cos 0er 2 + sin 9a x ) 
= 7 + 7-ct; 7= (j x , 7 Z ) 



where cr, are the Pauli matrices and A = (f a — /&)/2 = 
cos a. 

For the two-site Hubbard model in the sector of sin- 
glet states with N = 2 and S z = 0, (44 1 may be 
inverted to express = ^oM- Explicitly, we find 
*o = cos(a/2)<I> QQ — sin(a/2)<I> & b, where is the Slater 
determinant composed of the natural spin orbitals </>jj 
and 4>ii (i = a,b). The terms of the energy functional 
£[7] = T[ 7 ] + U[j\ + V[i\ are found to be 



Th] = 


(*o 


f 


*o) 


= -2£Asin0 


U[l] = 


(*o 


U 


*o) 




V[<y] = 


(*o 


V 




= V A cos0. 



sin 2 



(49) 



The electron-electron interaction energy functional U[y] 
agrees with t he ge neral exact result for 2-electron closed 
shell systems^ESI We may partition U [7] into the 
Hartree energy 

= C/(l + A 2 cos 2 
and the exchange-correlation energy 
E xc [l\ = U[i\-E H [j\ 

= -u (\ + Vi - a 2 



(50) 



U 
2" 



1 + A 2 + y/l -A 2 ) cos 



(51) 



In Sec. Ill Al the KS Hamiltonian was derived from the 
stationary principle for the energy. For the present model 
the KS Hamiltonian is a real 2x2 matrix. In the site 
basis its elements are h(ij) — (0\cihcj |0). This matrix 

may be expressed as h = h ■ a with 



B ( smao 

tlx = cos ao cos a sin ( 

4 \ sin a 

B sin ao(l + sin a) 2 . . , 

— : sm cos 1 

4 sin a cos a 



K = 



h, = 



B sinaofl + sin a) . 9 . 

; sin cos ( 

4 sin a cos a 



V 
2' 



(52) 



In these expressions the variable a represents the depen- 
dence on the occupation numbers through the definition 
a = cos^ 1 A = cos~ 1 ((/ a - / b )/2), a Q = tan _1 (C//4i) is 
the ground state value of a when V — 0, and represents 



the dependence on the natural orbitals, cf. (47 1. Let us 



verify (20) for the uniform case V — 0, for which the 
ground state 1-matrix has = n/2 and a = ao- At these 



values hr. 



(48) collapse in this case. 



0, which verifies the eigenvalue 
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2. Iteration of the KS equations 

We demonstrate here the iteration of the KS equa- 
tions following the straightforward algorithm described 
in Sec. |IIB| During the optimization of the orbitals the 
occupation numbers (i.e. a) are held fixed. Let us look 
more closely at each operation in the algorithm. In oper- 
ation (i) , t he K S Hamiltonian for step n + 1 is found by 



evaluating (52) at the 1-matrix 7 ("), i.e., at 8 = 8, 



operation (ii), we find the eigenstates Ui of M r 
we parametrize in the form (47 1 with 8 = 8 r< 



In 
which 
These 
~ 1 ' for 
equal 



eigenstates are taken as the natural orbitals 

step n + 1. This implies setting each of the cj)\ 
to one of the Ui . In the present case, the natural orbitals 
are chosen such that 8 n+ i is as close as possible to 8 n . In 
operation (iii), is constructed from the %™ by 

( 46 1 . We may now condense these three operations into 
a discrete iteration map on 8, i.e., a map 8 n — > 8 n+ \. It 
is defined by 



cos6>„+i = sgn(A - A gs ) 



(53) 



A aa is 



for < 8 n < 7T and A > (0 < a < tt/2). In (|53j, , 
the ground state value of A. An example of the iteration 
map for t = 1, U = 1, V = 0, and A = A gs — 0.02 is shown 
in Fig. [2] The solid [black] and dashed [red] curves are 
the left and right hand sides of ( 53 ) . The intersections of 
the two curves are fixed points of the iteration map. The 
ground state corresponds to the fixed point at 8 = tt/2. 
The iteration map may be represented graphically by al- 




FIG. 2: [color online] The iteration map (531 is shown for 
t = 1, U = 1, V = 0, and A = A qe - 0.02. The solid [black] 
curve is the left hand side of Q53D . The dashed [red] curve 
is the right hand side. The dotted [blue] curve demonstrates 
the first two iterations. The iteration map does not converge 
to the ground state fixed point 8 = 7r/2. 

ternately drawing vertical lines from the solid curve to 
the dashed curve and horizontal lines from the dashed 



curve to the solid curve. The dotted [blue] curve shows 
an example of the first two iterations beginning from an 
initial guess 8 . The next iterations 8\ and 82 move far- 
ther away from the ground state, and the map does not 
converge to the ground state fixed point 9 = tt/2. 

The iteration map is nonlinear and may exhibit quite 
complex behavior. The linearization of the map at a 
fixed point tells us whether the fixed point is stable or 
unstable. As an example, let us consider the uniform 
case V = 0, for which the ground state fixed point is 
8 = tt/2. Linearization of (53 1 in terms of the variable 
m = (ni — U2)/2 = A cos ( 



gives 



where 



m n+ i « sgn(A - A gs ) 



(1 + sin a) 2 
(cot ao — cot a) sin a cos a 



(54) 



(55) 



Suppose the occupation numbers 
ground state values, i.e., A = A 



<;.s 



are close to their 
+ 5 A where 6 A is 



a small displacement. The leading approximation for £ 



gives 



U 2 (U + B) 2 1 
TB3 SA' 



(56) 



For any nonzero values of t and U, there is a thresh- 
old d > such that for \SA\ < d, |£| > 1. Therefore 
the ground state is an unstable fixed point 



II B 



In Sec 

the divergence of the iteration map was connected to the 
divergence of the static KS response function. Let us 



verify (26) explicitly for the present case. As seen in 
(54 1, the linearized iteration map affects only the diag- 
onal elements of the 1-matrix, i.e., the density, which is 
described by the variable m. Therefore, the relevant re- 
sponse functions arc the density-density response for the 
KS system 



EE 



fi fj 



(4>i\m\ct>j)((t>j\m\(l>i) 



i jjti 
TU 3 1 

'Ipsa 



(57) 



and the density-density response for the interacting sys- 
tem 



X = 



E 

k 



^0 |m|*fc) |m|* 



0; 



Eq - E k 



+ c.c. 



'B(B + U)' 



(58) 



For the two-site Hubbard model, these response functions 
are just constants. The KS response has a functional de- 
pendence on the 1-matrix. It diverges as the ground state 
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is approached, i.e., in the limit 5A — > 0. The linearized 



iteration map ( 26 1 is simply multiplication by a constant 



XsX 



U 2 (U + B) 2 1 



TB 3, 



5A 



(59) 



which agrees with the direct calculation ( 56 I 



Of course, in actual calculations it is necessary to have 
a convergent iteration scheme. One possibility for ob- 
taining convergence is the level shifting method™ whose 
application in 1MFT was discussed in Sec. |II C| In the 
level shifting method, one introduces artificial shifts of 
the KS eigenvalues in order to improve convergence. A 
shift of the KS eigenvalue ei by an amount Q is equiva- 
lent to adding a term (i\4>i)(4>i\ to the KS Hamiltonian, 
where 4>i is the orbital with eigenvalue . The KS system 
for the two-site Hubbard model has two orbitals. As the 
divergence of the iteration map is due to the degeneracy 
of the KS spectrum at the ground state, it seems sensi- 
ble to prevent degeneracy by introducing a separation 2£ 
between the levels. Thus, we add the following term to 
the KS Hamiltonian at iteration step n 

-(\<fia)(<Pa\ + (\<t>b)(4>b\ = -C(shi0„a- x + cosd n a z ) , 

(60) 

where <p a and <pb are evaluated at — n . An example of 
the effect of level shifting is shown in Fig. [3j Convergence 




FIG. 3: [color online] The iteration map for t = 1, U = 5, 
V = —2.5, and A = A gs — 0.1. The solid [black] curve is 
the left hand side of (531. The (dotted [blue], dash-dotted 
[green] , dashed [red]) curve is the right hand side with level 
shift £ = (0, 3, 6). The threshold level shift for convergence is 
£ c ~ 4.07, which may be calculated with (1351). 



no physical meaning when the KS system is not self- 
consistent. The KS energy is shown in Fig. [4] as a func- 
tion of the deviation Sj = (5^ x ,Sj z ) of the 1-matrix 
(48 1 from the ground state 1-matrix. It is immediately 



Energy 




FIG. 4: [color online] The KS energy for t = 1, U = 3.5 
and V = is shown as a function of the deviation (S^/ X ,S"/ Z ) 
from the ground state. The two surfaces represent the two 
branches of the KS energy. The space curves show the energy 
as a function of 9 for fixed occupation numbers. Optimization 
of the orbitals corresponds to moving along one of these curves 
to the stationary point. The ground state is the point of conic 
intersection at the origin. 

seen that the KS energy is not stationary at the ground 
state 1-matrix, which is a cusp point where the energy 
E s changes linearly with respect to £7. The KS energy 
is multivalued due to the choice implied in occupying 
the KS levels when the system is not self-consistent (see 
Sec. II D ). The space curve in Fig. [1] shows the energy as a 
function of 9 for fixed occupation numbers, i.e., for fixed 
A. The KS response is proportional to the inverse sepa- 
ration between the two branches of the space curve. The 
separation vanishes as the curve approaches the conic 
point, which is the origin of the divergent KS response. 



C. Solution by DFT 



is achieved when £ exceeds a threshold, which may be cal- 
culated from the convergence criterion ( 35 1 . The dashed 



The two-site Hubbard model with the local external 



[red] curve in Fig. [3] shows the iteration map with a level 
shift value greater than the threshold. For the two-site 
Hubbard model, the criterion for convergence can be vi- 
sualized graphically as the condition that the magnitude 
of the slope of the level shifted curve be less than the 
slope of the solid [black] curve at the fixed point. 

At each iteration step the KS system has an "instan- 
taneous" energy E s [-f] = tr(h[j]j), which has, of course, 



potential chosen in (36) may be treated also with DFT. 



It is interesting to compare the DFT-KS scheme with 
the 1MFT-KS scheme, especially with regard to their 
convergence behavior. The variational energy functional 
and KS Hamiltonian may be constructed explicitly. An 
interesting result of the investigation is that the straight- 
forward iteration map is divergent when U > 1.307i (for 
V = 0). We derive a general condition for the conver- 
gence of the DFT-KS equations. 
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1. Energy functional 
The HK energy functional for a lattice is 

E[n,v]=^2v(t)ni + F[n], (61) 

i 

where v(i) is the external potential at site i and F[n] is a 
universal functional of the density (here, site occupancy) 
defined as 



F[n] = (V [n]\f + U\V [n}), 



(62) 



where T is the kinetic energy operator and U is electron- 
electron interaction. In the following treatment of the 
two-site Hubbard model, we depart from standard prac- 
tice by enforcing the normalization condition ni+ri2 = N 
explicitly (i.e., through the parametrization) , rather than 
with a Lagrange multiplier. Thus, we take as basic vari- 
able the single parameter m = {n\ — ri2)/2 that uniquely 
specifies the density (site occupancy). Similarly, the 
external potential is specified by the single parameter 
V = v(l) — v(2). The functional F[n] in (62 1 is then just 
a function F(m), which may be constructed explicitly as 
follows: i) a map m — > Wo is defined as the composition 
of the maps m — > V and V — > "Jo and ii) the result- 
ing function ^f (m) is used to evaluate (62 1. An explicit 
expression for the map m — ► V can be found from the in- 
verse of ( 43 ) . The second map v — > Wo was given in ( 40 1 . 



The composition of these two maps gives the ground state 
as a function of m, i.e., W (to), with which the universal 



functional (62 1 may be evaluated. 



2. KS Hamiltonian 

Following standard practice, the KS Hamiltonian takes 
over, unchanged, the kinetic energy operator from the 
many-body Hamiltonian. Thus, we consider the KS 
Hamiltonian 



h = t 



(63) 



where i = —t(c\c2 + c\c\) is the kinetic energy operator 
and v s (i) is the KS potential at site i defined by 



dW 

dm ' 



(64) 



where W[n] = E[n,v\ — T s [n] contains the Hartree and 
exchange-correlation energy as well as the external poten- 
tial energy, and T s [n] is the kinetic energy of the KS sys- 
tem. We do not separate these contributions explicitly. 
The KS potential is spin independent because the ground 
state density is spin unpolarized. Also, it is determined 
only to within an arbitrary additive constant, which we 
choose such that v s (l) + v s (2) = 0. In the site basis, the 
KS Hamiltonian is a 2 x 2 matrix which may be expressed 



as h — —ta x + (V s /2)a z , where = v s (l) — v s (2). The 
kinetic energy of the KS system is evaluated as 



T. = 



occ 



fi(4>i\t\(f>i) 



^{<Pa\t\4>a) 

-2tsin0, 



(65) 



where 4> a is the lowest energy eigenstate of (63 1 and is 



twice occupied (once by a spin up electron and once by 
a spin down electron.) It is parametrized as in (47) with 
tan# : 



-2t/V a . The density of the KS system is 



^2 fi(^i\m\<j>i) 



= cos 9. 



(66) 



Thus, from ( 65 1 and ( 66 ) the kinetic energy T s is a known 



function of m s . From (611, (64) and (65), the KS poten- 
tial is 



Vs 



dW 
dm 

d_ 
dm 

V + f(m s ) 



(E(m,V)-T s (m)) 
dm 



(67) 



where V — v(l) — v{2) is the given external potential and 



f(m s ) 



dF 
dm 



(68) 



Eq. 67 is simply the familiar expression v s (r) — v(r) + 
vh(t) + v xc (r) with a different partitioning of the terms. 
It is seen that the terms f — dT s / dm together correspond 
to the Hartree and exchange-correlation potentials. 



3. Iteration of the KS equations 

Let us investigate the iteration of the KS equations 
in the present context. The conventional iteration map 
consists of the following steps: i) the KS potential for 
step n + 1 is determined from the density of step n using 



(|67|, i.e., V} n+ ' = V s (mi n) ), ii) the eigenstates of h {n+1 ^ 
are found, and iii) the density of step n + 1 is calculated 
with (66). 



Consider step (i) in more detail. The KS potential is 
obtained from (luTl, 



=V + /(mW) 



dTs 
dm 



(69) 



The right hand side may be expressed differently by using 
the stationary conditions for the energy functional E[n, v] 
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and the KS energy E s = T s + v s (i)n 
condition dE/dra — applied to (61 1, gives / 



The stationary 
-V(m), 

where V(m) is the external potential such that the inter- 
acting system has ground state density m. Similarly, the 
stationary condition applied to E s gives dT s /dm = — V s . 
Substituting these relations in ( 69 1 yields 



1.0-1 



v (n+i) —y_ v( m W) + V s {m^). 



(70) 



At self-consistency the V s terms cancel, and we obtain the 
expected result V = V{m gs ), where m gs is the ground 
state density. For the present model, the ground state 
density could be found by solving V — V(m) as V(m) is 
known exactly from (43 ). However, in general the ground 
state must be found by iteration. Eq. [70] implies an it- 
eration map for the density, i.e. 



because V s 



(n+l) 



determines m 



(n+ 



a map m 
i) 



(n) 



(n+l) 



From ( 66 1 and the 



definition tan# = —2t/V s , we find the relationship 



V s 



-2t 



(71) 



The density may be iterated until self-consistency is 
reached. However, we encounter a technical difficulty 
for the present mode l. In order to express explicitly the 
term V{m s n) ) in (70 I, we must invert (W3| , which involves 
solving a cubic equation. As the solutions are rather 
unwieldy, we take here a different approach. We iterate 
instead the external potential V(m). It may seem odd 
to iterate the external potential, which is given in the 
statement of the problem. Nevertheless, the iteration 
map for V provides an "image" of the iteration map for 
m s , by virtue of the HK theorem. Such an approach 
allows us to investigate certain features of the iteration 
map, in particular its convergence behavior. In order to 
express ( 70 1 as an iteration map for V, we need to express 
V s as a function of V. In other words, we find the value of 
V s such that the KS system has density m s — m, where 
77i is the density of the interacting system with V. The 
composition of (71 1 and d43| yields the desired function 



V S {V) = V s {m{V)). 



(72) 



Using (72 1 in (70 1, we obtain the iteration map for the 

(73) 



external potential 

V s (V (n+1) ) =V- V (n) + V s {V {n) ), 



which is expressed in implicit form. 

Examples of the iteration map for a uniform system 
(V = 0) are shown in Figs. [5] and 6] where the left and 
right hand sides of (73 1 are plotted. Suppose an initial 
value V^ ) ^ is chosen. For a system with V = 0, the 
ground state has uniform density (m = 0), but the initial 
density mi°^ associated with is not uniform. Upon 
iteration, we expect the KS system to relax to a uniform 
density, i.e., we expect the KS potential to be such as to 
push the system closer to uniform occupancy in the next 
iteration. The solid [black] curves in Figs. [5] and [6] rep- 
resent the left hand side of (73 1, while the dashed [red] 




-1.0- 1 



FIG. 5: The iteration map for the external potential V is 
shown for a weakly interacting system with U — t. The left 
and right hand sides of ( 73 1 are shown as solid [black] and 
dashed [red] curves, respectively. 




-1.0— 1 



FIG. 6: The iteration map for the external potential V is 
shown for a strongly interacting system with U = 4t. The 
left and right hand sides of ( 73 1 are shown as solid [black] 



and dashed [red] curves, respectively. 



curves represent the right hand side. The iteration map 
may be demonstrated graphically by alternately drawing 
vertical lines from the solid curve to the dashed curve and 
horizontal lines from the dashed curve to the solid curve. 
The map displays "charge oscillation." The ground state 
is a stable fixed point if the magnitude of the slope of 
the dashed curve at the origin is less than the slope of 
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the solid curve at the origin. For weakly interacting sys- 
tems the iteration map is convergent, while for strongly 
interacting systems it is nonconvergent . The threshold 
for convergence is U « 1.307£. 

4- Linearization of the KS equations 



is the exchange-correlation kernel. The necessary and 
sufficient condition for convergence of the KS equations 
is that all eigenvalues of the operator 

Xs{vc + Lc) (80) 
have modulus less than 1. 



The nature of the fixed point and the origin of diver- 
ent behavior arc revealed by linearization of the iteration 
map. We linearize the map by expanding both sides of 
( 70 1 with respect to Sm s = m s — m 



gs , where m gs is the 



ground state density. We find 



(W" +1 > 



X, 



X 



8m { T 



(74) 



where Xs and x are the density-density response func- 
tions defined in ( 57 1 and ( 58 1 . The threshold for conver- 
gent behavior is 



XsX 



< 1, 



(75) 



or equivalently, XsX 1 < 2. Note the change from 1 for 
1MFT to 2 for DFT on the right hand side, cf. ((26 1. 
Consider the case V — 0, which has uniform density in 
the ground state (m gs = 0). Using the (57 1 and (58 1 in 
( 75 1 , gives the threshold condition 



cos(?r/4 - a /2) = 4 (sin(7r/4 - a /2)Y 



(76) 



The threshold is cos(7r/4 — ao/2) ~ 0.8095, which corre- 
sponds to U w 1.307£. Let us consider the limit U — > oo. 
The leading behavior of the KS response is independent 
of U, 



1 



(77) 



while the response of the interacting system vanishes as 



X 



III 
4[73' 



(78) 



For sufficiently large U, the threshold ( 75 ) is crossed and 
divergent behavior results. In DFT, as also in 1MFT, 
the source of divergent behavior is a KS response that is 
too large in relation to the exact response. In 1MFT the 
imbalance results from a divergent KS response, whereas 
in DFT the KS response generally remains finite but the 
response of the interacting system becomes too small as 
U increases. 

In standard DFT (with continuous n(r)), the analog 
of the linearized iteration map ( 74 ) may be written 



n {n+1 \r) » J dr'dr"x s (ry)(v c (r , ,r") + f xc (r , y')) 

x»< n V)> ( 79 ) 

where n^(r) is the density of iteration step n, v c is the 
kernel of the Coulomb interaction, and f xc = Sv xc /5n 



IV. CONCLUSIONS 

The status of the KS system in 1MFT has been uncer- 
tain. Although Gilbert derived effective single-particle 
equations from the stationary conditions for the energy 
functional, the degeneracy of essentially all of th e result- 
ing orbitals was thought to be paradoxical! 2 * 13 * 1 ^ We have 
presented an alternative derivation of the KS equations 
in which the degeneracy is lifted by constraining the oc- 
cupation numbers. Such a KS scheme is well-behaved in 
the neighborhood of the ground state occupation num- 
bers. Therefore, the correct natural orbitals are obtained 
in the limit that the ground state is approached. We have 
constructed explicitly the IMFT-KS system for a simple 
two-site Hubbard model. While we find no paradoxi- 
cal results, the KS system has many striking features, in 
particular the collapse of eigenvalues at the ground state. 
Although the KS eigenvalues do not have a physical in- 
terpretation as in DFT, the orbitals, which are called 
natural orbitals, play an important role in the context 
of configuration interaction, i.e., the expansion of the 
full wavefunction as a sum of Slater determinants. 1 ' This 
may be important in the search for approximate energy 
functionals. 

Beyond the question of the existence of the KS system 
in 1MFT, there is the issue of its practicality. The KS 
system has been extremely useful in DFT calculations. 
Due to the implicit 1-matrix dependence of the single- 
particle potential, the KS equations are nonlinear. Such 
equations are generally solved by iteration. As in DFT, 
there is a "straightforward" procedure for iteration. In 
contrast to DFT, the "straightforward" procedure is al- 
ways divergent, in the sense that the ground state is an 
unstable fixed point. We have demonstrated the insta- 
bility of the ground state by linearization of the itera- 
tion map. The source of the instability is the divergence 
of the KS static response function at the ground state, 
which in turn, is due to the degeneracy of the KS spec- 
trum. Degeneracy-driven instability is reminiscent of the 
Jahn-Teller effect, and the connection is strengthened if 
we regard the implicit 1-matrix dependence of the KS 
Hamiltonian as analogous to the parametric dependence 
of the Born-Oppenheimer Hamiltonian on nuclear coordi- 
nates. In both cases, the energy spectrum splits linearly 
with respect to displacement from the degeneracy point. 
Thus, the energy may always be lowered by displacement. 
For the IMFT-KS system, this means that the KS energy 
trijvy) may always be lowered by displacement from the 
ground state, leading to an instability of the iteration 
procedure. However, this is a fictitious energy and the 
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HK energy functional E v is of course always minimum at 
the ground state. 
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APPENDIX A: GROUND STATE NOT 
DETERMINED BY THE DENSITY 

We give here a simple example which shows that the 
density alone does not always uniquely determine the 
ground state wavefunction if the external potential is 
nonlocal. Our example is the two-site Hubbard model, 
which was so lved for the case of a local external potential 



Sec. 



Ill A 



The Hamiltonian H = T + V + U is given 
in ( 36 1 . In such a lattice model, the hopping parameters 
tjk are real numbers that represent the kinetic energy. A 
"magnetic field" can be introduced by giving tjh a phase, 
i.e., by the transformation tjk — > £jfcexp(t2 n -v4.(n)), 
where A(n) is the "vector potential" at site n and the sum 
runs over a string of sites from site j to site k. For the 
two-site model this is just the transformation t — > te lT . 
We see that this magnetic field appears in the Hamil- 
tonian in exactly the same manner as a nonlocal exter- 
nal potential, such as V\2c\c2 + h.c, because it modifies 
the nonlocal hopping terms. We can generate the above 
phase transformation by the U(l) rotations c\ — > cie lT / 2 



and C2 — > Ci& 



-ir/2 



The cigenstates of the transformed 



Hamiltonian are readily generated from the eigenstates 
of the original Hamiltonian by applying the same trans- 
formation. For example, without the magnetic field, the 
ground state to first order in small V is 



<$ 2 |t>|$ ) 



$9 



(Al) 



where $j are given in (37 1. When the magnetic field is 



turned on, the $j change, e.g. 




(A2) 



in the site basis 



2T 1| 



o 2 |^ 2 l |0) }. Accordingly, the ground state acquires 
a nontrivial dependence on the magnetic field (r- 
dependence). At the same time, the ground state 1- 
matrix is transformed to 



7 = 70 + A Q 



COS Uq 

sin Bqc 1 



smtfoe 
— cos 9q 



(A3) 



where Aq and 9q are the ground state values (for r = 0) 
of the variables defined in ( 47 ) and ( 48 ) . The density is 



given by the diagonal elements, which are unaffected by 
the transformation. Only the off-diagonal (nonlocal) ele- 
ments are sensitive to the magnetic field. Therefore, the 
1-matrix rather than the density is required to uniquely 
specify the ground stateP 
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The following theorenP^ implies that the functional E v [7] 
is stationary, i.e., SE V — with respect to all variations 
in the VR space. Theorem - A necessary condition for the 
differentiable functional J[y] to have an extremum (mini- 
mum) for y — yo is that its first variation vanish for y — yo, 
i.e., that 5J[h] — for y = yo and all admissible variations 
h. Thus, the functional E v [7] must be stationary because it 
is differentiable in the VR space and the quantum mechani- 
cal variation principle (Rayleigh-Ritz variational principle) 
implies that it is minimum at the ground state 1-matrix 
7 = "/ gs . The differentiability of E v [7] follows from the dif- 
ferentiability of the mapping 7 — > >I<o in the VR space, if 
the ground state \l/o is nondegenerate. 
If the DFT-KS system is degenerate, the occupation num- 
bers of the degenerate KS orbitals are not determined by 
the aufbau principle. In this case, the KS system adopts 
an ensemble state, and the occupation numbers of the de- 
generate orbitals are chosen such that the KS system is 
self-consist ent a nd reproduces the density of the interact- 
ing systempIEU 

We do not demonstrate this statement for \ as an operator 
on an infinite dimensional space. However, it is valid in a 
finite dimensional basis (see Ref. 1331 for a proof in DFT). 
For convenience we assume here that the system has no 
pinned states. 



